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made available via the internet — there are no plans to publish it. 

Abstract 

This is a introduction to the Jacobian transpose method, the 
pseudoinverse method, and the damped least squares methods for 
inverse kinematics (IK). The mathematical foundations of these 
methods are presented, with an analysis based on the singular value 
decomposition. 


1 Introduction 

A rigid multibody system consists of a set of rigid objects, called links, 
joined together by joints. Simple kinds of joints include revolute (rotational) 
and prismatic (translational) joints. It is also possible to work with more 
general types of joints, and thereby simulate non-rigid objects. Well-known 
applications of rigid multibodies include robotic arms as well as virtual 
skeletons for animation in computer graphics. 

To control the movement of a rigid multibody it is common to use inverse 
kinematics (IK). For IK, it is presumed that specified points, called “end 
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effectors,” on the links are are assigned “target positions.” To solve the 
IK problem, we must find settings for the joint angles so that the resulting 
configuration of the multibody places each end effector at its target position. 
More general formulations of IK allow also orientation goals, or directional 
goals. 

There are several methods for solving IK problems, coming originally 
from robotics applications. These include cyclic coordinate descent meth¬ 
ods [43], pseudoinverse methods [45], Jacobian transpose methods [5, 46], the 
Levenberg-Marquardt damped least squares methods [41, 34], quasi-Newton 
and conjugate gradient methods [43, 49, 15], and neural net and artificial 
intelligence methods [19, 27, 36, 38, 20, 22, 40, 16]. 

The present paper focuses on applications of IK in computer graphics and 
real-time animation. There has already been extensive use of IK in computer 
graphics [18, 26, 25, 44, 23, 2, 1, 17, 24, 29, 39, 21, 37, 12]: the most common 
applications are animating humans or creatures by specifying the positions, 
and possibly the orientations, of their hands, feet and head. Our interests lie 
particularly in using target positions for end effectors to animate an entire 
multibody, and in methods that are robust and behave well in wide range of 
situations. As part of the robustness, we want the end effectors to track the 
target positions and to do a reasonable job even when the target positions are 
in unreachable positions. In this paper, we consider only first order methods 
and consider the following generic application: we presume a multibody has 
multiple end effectors and multiple target positions, given in real-time in 
an online fashion, and want to update the multibody configuration so as to 
dynamically track the target positions with the end effectors. 

One might wonder why it is important to allow target positions to 
be unreachable. There are several reasons: First, it may be difficult to 
completely eliminate the possibility of unreachable positions and still get 
the desired motion. Second, if target positions are barely reachable and 
can be reached only with full extension of the links, then the situation is 
very similar to having unreachable targets. Unfortunately, the situation 
of target positions in unreachable positions is difficult to handle robustly. 
Many methods, such as the pseudoinverse or Jacobian transpose methods, 
will oscillate badly in this situation; however, (selectively) damped least 
squares methods can still perform well with unreachable target positions. 

The outline of the paper is as follows. We first introduce the mathemat¬ 
ical framework for the IK problem. We then discuss the Jacobian transpose 
method, the pseudoinverse method, the singular value decomposition, and 
the damped least squares (DLS) method. For an extension of the DLS 
methods to a method called selectively damped least squares (SDLS), see 
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Buss and Kim [7]. Nearly all the present paper is expository, but new 
aspects include the possibility of forming the Jacobian matrix with the 
target positions instead of the end effector positions. We attempt to explain 
the mathematical foundations clearly so as to elucidate the strengths and 
weaknesses of the various methods. 

For simplicity and to keep the paper short, we do not consider any 
aspects of joint limits or avoiding self-collisions; rather, we will only consider 
the “pure” IK problem without joint limits and without self-collisions. 

2 Preliminaries: forward kinematics and Jacobians 

A multibody is modeled with a set of links connected by joints. There 
are a variety of possible joint types. Perhaps the most common type is 
a rotational joint with its configuration described by a single scalar angle 
value. Other joint types include prismatic (he., translational, or sliding) 
joints, screw joints, etc. For simplicity, we will discuss only rotational joints, 
but the algorithms and theory all apply to arbitrary joints. The key point is 
that the configuration of a joint is a continuous function of one or more real 
scalars; for rotational joints, the scalar is the angle of the joint. 

The complete configuration of the multibody is specified by the scalars 
9\,6 n describing the joints’ configurations. We assume there are n joints 
and each Oj value is called a joint angle (but, as we just said, could more 
generally represent a value which is not an angle). Certain points on the links 
are identified as end effectors. If there are k end effectors, their positions 
are denoted si,... ,Sfc. Each end effector position Sj is a function of the 
joint angles. We write s for the column vector (si,S 2 ,..., Sf ; ) T ; this can be 
viewed as a column vector either with m = 3k many scalar entries or with 
k many entries from R 3 . 

The multibody will be controlled by specifying target positions for 
the end effectors. The target positions are also given by a vector 
t = (ti,... ,tfc) T , where t* is the target position for the ith end effector. 
We let e* = tj — Sj, the desired change in position of the ith end effector. 
We also let e = t — s . 

The joint angles are written as a column vector as 6 = (#i,..., 0 n ) T . 
The end effector positions are functions of the joint angles; this fact can be 
expressed as s = s(0 ), or, for i = 1..... fc, s, = s,(0). The IK problem is to 
find values for the 9j ’s so that 

tj = s i(0), for alii. (1) 
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Unfortunately, there may not always be a solution to (1), and there may not 
be a unique (best) solution. Even in well-behaved situations, there may be 
no closed form equation for the solution. 

We can, however, use iterative methods to approximate a good solution. 
For this, the functions sj are linearly approximated using the Jacobian 
matrix. The Jacobian matrix J is a function of the 9 values and is defined 
by 



Note that J can be viewed either as a kxn matrix whose entries are vectors 
from M 3 , or as m x n matrix with scalar entries (with rn = 3k). 

The basic equation for forward dynamics that describes the velocities 
of the end effectors can be written as follows (using dot notation for first 
derivatives): 

S=J(0)0. (2) 

The Jacobian leads to an iterative method for solving equation (1). 
Suppose we have current values for 9, s and t. From these, the Jacobian 
J = J{9) is computed. We then seek an update value A9 for the purpose 
of incrementing the joint angles 0 by A 6: 

9 := 9 +A9. (3) 

By (2), the change in end effector positions caused by this change in joint 
angles can be estimated as 

As & JA0. (4) 

The idea is that the A 9 value should chosen so that As is approximately 
equal to e, although it is also common to choose A 9 so that the approximate 
movement As in the end effectors (partially) matches the velocities of the 
target positions (see [45]). The update of the joint angles can be used in two 
modes: (i) Each simulation step performs a single update to the value of joint 
angles using equation (3), so that the end effector positions approximately 
follow the target positions, (ii) The joint angles are updated iteratively until 
a value of s is obtained that is sufficiently close to a solution. It is also 
possible to use a hybrid of (i) and (ii), that is, using a small number of 
repeated updates using (3) so as to more accurately track the end effector 
positions. 

The rest of this paper discusses strategies for choosing A 9 to update the 
joint angles. In light of (4), one approach is to solve the equation 


JA9 
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The entries in the Jacobian matrix are usually very easy to calculate. 
If the j th joint is a rotational joint with a single degree of freedom, the 
joint angle is a single scalar 9j . Let p j be the position of the joint, and 
let Vj be a unit vector pointing along the current axis of rotation for the 
joint. In this case, if angles are measured in radians with the direction 
of rotation given by the right rule and if the i th end effector is affected 
by the joint, then the corresponding entry in the Jacobian is 


ds i 
Wj 


vj x (sj - p|| 


If the ith end effector is not affected by the j th joint, then of course 
(9s j / 99j = 0. 

If the jth joint is translational, the entry in the Jacobian matrix is 
even easier to compute. Suppose the j th joint performs translation the 
direction of the unit vector Vj , so that the the joint “angle” measures 
distance moved in the direction Vj. Then if the ?'th end effector is 
affected by the jth joint, we have 


ds i 

Wj 


For more information, see Orin and Schrader [35] who discuss how 
to calculate the Jacobian matrix entries for different representations 
of joints and multibodies. The textbook [6, Ch. 12] also discusses a 
representation of rigid multibodies and how to calculate the Jacobian. 

Calculating the Jacobian 


for A 6 . In most cases, this equation cannot be solved uniquely. Indeed, the 
Jacobian J may not be square or invertible, and even if is invertible, just 
setting A 6 = J -1 e may work poorly if J is nearly singular. 

An alternate Jacobian. An alternate method for defining the Jacobian 
matrix is to let 



where the partial derivative is calculated using the formula for ( 8si/89j ) 
with t i substituted for s,;. The meaning of 8t l j89 ] is that the target 
position is thought of as being attached to the same link as the i th end 
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effector. The intuition is that with this formulation of the Jacobian, we are 
trying to move the target positions towards the end effectors, rather than 
the end effectors towards the target positions. 

The alternate Jacobian may be used in place of the usual Jacobian in 
any of the algorithms discussed below. Our experience has been that this 
alternate can improve on the usual Jacobian in terms of reducing oscillation 
or overshoot when target positions are too far away to be reached by the 
end effectors. However, the drawback is that in some configurations, the 
alternative Jacobian can lead to “jerky” behavior. This is particularly true 
for rotational joints when the multibody’s links are folded back on each other 
trying to reach a close target position. 

Setting target positions closer. A recurring problem in tracking target 
positions, is that when the target positions are too distant, the multibody’s 
arms stretch out to try to reach the target position. Once the multibody is 
extended in this way, it usually is near a singularity (that is, the Jacobian is 
very sensitive to small changes in joint angles), and the multibody will often 
shake or jitter, attempting unsuccessfully to reach the distant target. These 
effects can be reduced with DLS and SDLS algorithms, but are difficult to 
remove completely. 

One technique to reduce this problem is to move the target positions in 
closer to the end effector positions. For this, we change the definition of e; 
instead of merely setting e = t — s , each component e,; in the vector e has 
its length clamped to a specified maximum value. That is, we define 

e* = ClampMag(tj — Sj, T> max ), 

where 

r w if i|w|| < d 
ClampMag(w, d) = < 

( otherwise 

Here ||w|| represents the usual Euclidean norm of w. The value D max is 
an upper bound on how far we attempt to move an end effector in a single 
update step. 

For damped least squares, clamping the magnitudes of e* in this way 
can reduce oscillation when target positions are out of reach. This has 
the advantage of allowing the use of a smaller damping constant; the 
smaller damping constant allows significantly quicker convergence to target 
positions.* When the end effectors are tracking continuously moving target 

*Also for SDLS, [7] have found that clamping the magnitudes of e, in this way can 
effectively reduce oscillation when target positions are out of reach. 
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positions, the An ax distance should be at least several times larger than an 
end effector moves in a single update step. In our experience, setting D max 
to be approximately half the length of a typical link works well. 

For target positions that may jump discontinuously, we have used 
separate maximum values D maxi for each i . After a discontinuous movement 
of the target positions (or when beginning a simulation of a continuously 
moving target), we initially set D maxi to infinity. After the first simulation 
step, we let di be the amount by which the previous simulation step 
moved the ith end effector closer to its target position. Then, we let 
T>max, i = di + An ax , and use -D ma x,i to clamp the magnitude of . 

3 The Jacobian transpose method 

The Jacobian transpose method was first used for inverse kinematics by 
[5, 46]. The basic idea is very simple: use the transpose of J instead of the 
inverse of J. That is, we set Ad equal to 


a J 1 e, 


for some appropriate scalar a. Now, of course, the transpose of the Jacobian 
is not the same as the inverse; however, it is possible to justify the use of the 
transpose in terms of virtual forces. More generally, it can be shown that 
the following theorem holds [5, 46]. 

Theorem 1 For all J and e, (JJ T e,e) > 0. 

Proof The proof is trivial: (JJ T e, e) = (J T e,J T e) = ||J T e|| 2 >0. □ 


The approximation (4) implies that, for sufficiently small a > 0, 
updating the angles by equation (3) using A 6 = aJ T e will change the 
end effector positions by approximately aJJ T e. By Theorem 1, this has the 
effect of reducing the magnitude of the error vector e if a is small enough. 

It remains to decide how to choose the value of a. One reasonable way 
to try to minimize the new value of the error vector e after the update. 
For this, we assume that the change in end effector position will be exactly 
aJJ T e, and choose a so as to make this value as close as possible to e. 
This gives 


a 


(e, J J T e) 
(JJ T e, JJ T e )' 
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4 The pseudoinverse method 

The pseudoinverse method sets the value A 6 equal to 

A 6 = J% (6) 

where the n x m matrix jt is the pseudoinverse of J, also called the 
Moore-Penrose inverse of J. It is defined for all matrices J, even ones 
which are not square or not of full row rank. The pseudoinverse gives 
the best possible solution to the equation JAG = e in the sense of least 
squares. In particular, the pseudoinverse has the following nice properties. 
Let A G be defined by equation (6). First, suppose e is in the range (i.e., the 
column span) of J. In this case, JAG = e; furthermore, AO is the unique 
vector of smallest magnitude satisfying JAG = e. Second, suppose that e 
is not in the range of J. In this case, JAG = e is impossible. However, 
AG has the property that it minimizes the magnitude of the difference 
JAG — e. Furthermore, AG is the unique vector of smallest magnitude 
which minimizes 11 JAG — e| | , or equivalently, which minimizes | \JAG — e| | 2 . 

The pseudoinverse tends to have stability problems in the neighborhoods 
of singularities. At a singularity, the Jacobian matrix no longer has full row 
rank, corresponding to the fact that there is a direction of movement of the 
end effectors which is not achievable. If the configuration is exactly at a 
singularity, then the pseudoinverse method will not attempt to move in an 
impossible direction, and the pseudoinverse will be well-behaved. However, 
if the configuration is close to a singularity, then the pseudoinverse method 
will lead to very large changes in joint angles, even for small movements in 
the target position. In practice, roundoff errors mean that true singularities 
are rarely reached and instead singularity have to be detected by checking 
values for being near-zero. 

The pseudoinverse has the further property that the matrix (I — jt J) 
performs a projection onto the nullspace of J. Therefore, for all vectors <p, 
J(I — J)ip = 0. This means that we can set AG by 

AG = J^e + (/ — J^J)<p (7) 

for any vector p and still obtain a value for AG which minimizes the value 
JAG — e. This nullspace method was first exploited Liegeois [28], who used 
it to avoid joint limits. By suitably choosing ip, one can try to achieve 
secondary goals in addition to having the end effectors track the target 
positions. For instance, ip might be chosen to try to return the joint angles 
back to rest positions [18]: this can help avoid singular configurations. 







A number of authors (see [4]) have used the nullspace method to help 
avoid singular configurations by maximizing Yoshikawa’s manipulability 
measure [48, 47]. Maciejewski and Klein [30] used the nullspace method 
for obstacle avoidance. A more sophisticated nullspace method, called the 
extended Jacobian method , was proposed by Baillieul [4]: in the extended 
Jacobian method a local minimum value of a function is tracked as a 
secondary objective. The nullspace method has also been used to assign 
different priorities to different tasks (see [10, 3]). 

An algorithm for the pseudoinverse method can be derived as follows: 
From equation (5), we get the normal equation 

J T J AO = J T e. 

Then we let z = J T e and solve the equation 

(J t J)A9 = z. (8) 

Now it can be shown that z is always in the range of J T J , hence equation (8) 
always has a solution. In principle, row operations can be used to find the 
solution to (8) with minimum magnitude; however, in the neighborhood of 
singularities, the algorithm is inherently numerically unstable. 

When J has full row rank, then JJ T is guaranteed to be invertible. 
In this case, the minimum magnitude solution AO to equation (8) can be 
expressed as 

AO = J T (JJ T )- 1 e. (9) 

To prove this, note that if AO satisfies (9), then AO is in the row span of J 
and JAO = e. Equation (9) cannot be used if J does not have full row 
rank. A general formula for the pseudoinverse for J not of full row rank can 
be found in [6]. 

The pseudoinverse method is widely discussed in the literature but 
it often performs poorly because of instability near singularities. The 
(selectively) damped least squares methods have much superior performance. 

5 Damped least squares 

The damped least squares method avoids many of the pseudoinverse 
method’s problems with singularities and can give a numerically stable 
method of selecting AO. It is also called the Levenberg-Marquardt method 
and was first used for inverse kinematics by Wampler [41] and Nakamura 
and Hanafusa [34]. 




The damped least squares method can be theoretically justified as 
follows (see [42]). Rather than just finding the minimum vector A 6 that 
gives a best solution to equation (5), we find the value of Ad that minimizes 
the quantity 

\\JA0 — e|| 2 + A 2 ||A0|| 2 , 

where A € K is a non-zero damping constant. This is equivalent to 
minimizing the quantity 


( J ^ Ad - 

( s ) 

yxi) 

w 


The corresponding normal equation is 



This can be equivalently rewritten as 

(J T J + X 2 I)A6 = J T e. 

It can be shown (by the methods of section 6 below) that J T J + A 2 / is 
non-singular. Thus, the damped least squares solution is equal to 

A 6 = (J T J + A 2 !)" 1 /^. (10) 

Now J T J is an n x n matrix, where n is the number of degrees of freedom. 
It is easy to show that ( J T J + A 2 I)~ l J T = J T (JJ T + A 2 /) -1 . Thus, 

A 6 = J T (JJ T + \ 2 I)~ 1 e. (11) 

The advantage of equation (11) over (10) is that the matrix being inverted is 
only m x m where m = 3k is the dimension of the space of target positions, 
and m is often much less than n. 

Additionally, (11) can be computed without needing to carry out the 
matrix inversion, instead row operations can find f such that (JJ T +A 2 /)f = 
e and then J T f is the solution. 

The damping constant depends on the details of the multibody and 
the target positions and must be chosen carefully to make equation (11) 
numerically stable. The damping constant should large enough so that the 
solutions for Ad are well-behaved near singularities, but if it is chosen too 
large, then the convergence rate is too slow. There have been a number of 
methods proposed for selecting damping constants dynamically based on the 
configuration of the articulated multibody [34, 14, 15, 31, 11, 13, 8, 9, 33, 32]. 
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6 Singular value decomposition 


The singular value decomposition (SVD) provides a powerful method for 
analyzing the pseudoinverse and the damped least squares methods. In 
addition, we shall use the SVD to design a selectively damped least squares 
method in [7]. Let J be the Jacobian matrix. A singular value decomposition 
of J consists of expressing J in the form 

J = udv t , 

where U and V are orthogonal matrices and D is diagonal. If J is rri x n, 
then U is m x m, D is mx n, and V is n x n. The only non-zero entries 
in the matrix D are the values o % = d % , % along the diagonal. We henceforth 
assume m < n. Without loss of generality, a\ > <72 > • • • > a m > 0. 

Note that the values cq may be zero. In fact, the rank of J is equal 
to the largest value r such that a r ^ 0. For i > r, a t = 0. We use u, 
and Vi to denote the zth columns of U and V. The orthogonality of U 
and V implies that the columns of U (resp., of V ) form an orthonormal 
basis for M rn (resp., for R n ). The vectors v r+ i,..., v n are an orthonormal 
basis for the nullspace of J. The singular value decomposition of J always 
exists, and it implies that J can be written in the form 

Y^aiUivf. (12) 

The transpose, D T , of D is the n x m diagonal matrix with diagonal 
entries cq = d lyl . The product DD T is the m x m matrix with diagonal 
entries . The pseudoinverse, = (d\ -), of D is the n x m diagonal 
matrix with diagonal entries 

1 /dij if di.i d 0 
0 ’ if di,i = 0. 

The pseudoinverse of J is equal to 

Jt = VD ] U t . 



J = = 


Thus, 

jt = E^rV- 


(13) 
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The damped least squares method is also easy to understand with the 
SVD. The matrix JJ T + A 2 1 is equal to 

JJ t + A 2 I = (UDV t )(VD t U t ) + A 2 / = U(DD T + \ 2 I)U T . 

The matrix DD T + A 2 1 is the diagonal matrix with diagonal entries erf + A 2 . 
Clearly, DD T + A 2 1 is non-singular, and its inverse is the m x m diagonal 
matrix with non-zero entries (erf + A 2 ) -1 . Then, 

J T (JJ T + A 2 /) -1 = {VD T (DD T + \ 2 I)~ 1 U T = VEU t , 

where E is the n x m diagonal matrix with diagonal entries equal to 

- - ^ 

Thus, the damped least squares solution can be expressed in the form 

j t (jj t + a 2 /)- 1 = (M) 

Comparison of equations (13) and (14) makes clear the relationship 
between the pseudoinverse and damped least squares methods. In both 
cases, J is “inverted” by an expression T,y t uJ . For pseudoinverses, the 
value Ti is just a^ 1 (setting 0 -1 = 0); whereas for damped least squares, 
Tj = <7i/ (erf + A 2 ). The pseudoinverse method is unstable as a, approaches 
zero; in fact, it is exactly at singularities that cq’s are equal to zero. 

For values of a t which are large compared to A, the damped least squares 
method is not very different from the pseudoinverse since for large cq, 
cq/(crf + A 2 ) « 1 / a,. But, when o % is of the same order of magnitude as A 
or smaller, then the values cr _1 and (Ti/(af + A 2 ) diverge. Indeed, for any 
A > 0, cq/(crf + A 2 ) —> 0 as <jj —► 0. Thus, the damped least squares method 
tends to act similarly to the pseudoinverse method away from singularities 
and effectively smooths out the performance of pseudoinverse method in the 
neighborhood of singularities. 

7 Selectively damped least squares 

For the material that used to be in this section, see Buss and Kim [7]. 


12 





Figure 1: The Y and double-Y shapes. The end effectors are at the ends of 
the branches; the red balls indicate the target positions. 

8 Experimental results and recommendations 

For the rest of the material that used to be in this section, plus the results of 
additional experiments, see Buss and Kim [7]. 

To compare the IK algorithms, we implemented the “Y”-shaped and 
“double-Y” shaped multibodies pictured in figure 1. The first has seven 
links with two end effectors and the latter has 16 links with 4 end effectors. 
We let the target positions (the red balls in the figures) move in sinusoidally 
varying curves in and out of reach of the multibodies. The target positions 
moved in small increments (just large enough to still look visually smooth), 
and in each time step we updated the joint angles once.t Since joint angles 
were updated only once per time step, the end effectors tracked the target 
positions only approximately, even when the target positions were within 
reach. We visually inspected the simulations for oscillations and tracking 
abilities. We also measured the accuracy of the tracking over a period of 
hundreds of simulation steps. 

The Jacobian transpose had the advantage of being fast, but of poor 
quality. Its quality was poor for the Y shape and extremely poor for the 
double-Y shape. However, in other simultations, we have seen the Jacobian 
transpose method work well for a system with a single end effector. 

The pseudoinverse method worked very poorly whenever the target 
positions were out of reach, and we do not recommend its use unless joint 

^All our software, including source code, is available from the web page 
http://math.ucsd.edu/~sbuss/ResearchWeb/ikmethods. Short movie clips of the 
Jacobian transpose, the pure pseudoinverse method, the DLS and the SDLS methods are 
also available there. 
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angles are severely clamped with ClampMaxAbs. 

The damped least squares method worked substantially better than the 
Jacobian transpose method, although it is somewhat slower. We attempted 
to set the damping constant A so as to minimize the average error of the end 
effectors’s positions, but at the point where the error was minimized, there 
was a lot of oscillation and shaking. Thus, we had to raise the damping 
constant until unwanted oscillation became very rare (but at the cost of 
accuracy in tracking the target positions). 

We also implemented a version of the DLS method which uses the 
ClampMag method to clamp the components of the e vector: this method 
is called DLS'. The advantage of the DLS' method is that the clamping of e 
reduces oscillation and shaking, and thus a lower damping constant can be 
used. The lower damping constant allows the multibody to more agressively 
move towards the target positions. 

The runtimes for two different methods are described in the table below. 
Runtimes are in microseconds and were measured with custom C++ code on 
a 2.8GHz Pentium. The DLS' runtime is not reported, but is very close to 
that of DLS. For the Y-shape, the Jacobian matrix is 6x7, for the double-Y, 
it is 12 x 16. 


Shape 

Jacobian 

Transpose 

DLS 

Y 

1.1 fiS 

2.2 /j,s 

Double-Y 

6.5 /as 

18.5 /u.s 


We conclude with some recommendations. First, the Jacobian transpose 
performed poorly in our tests, but we have seen it work well in situations 
where there is a single end effector. For these applications, the Jacobian 
transpose is fast and easy to implement. For multiple end effectors, the DLS 
or DLS' methods can be used. For controlled situations where a damping 
constant can be set ahead of time, the DLS' method gives good performance 
and relatively easy implementation. For recommendations relating to the 
use of selectively damped least squares, see [7]. 

A final recommendation that applies to any method is that it is almost 
always a good idea to clamp the maximum angle change in a single update 
to avoid bad behavior from unwanted large instantaneous changes in angles. 
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